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^ Contributions  to  solving  initial  boundary  value  problems  for  partial 
differential  equations  have  been  made  by  applying  finite-difference  methods 
to  solve  seismic  wave  propagation  problems.  Very  little  has  been  done  in 
the  area  of  underwater  acoustic  wave  propagation  problems,  although  a set 
of  properly  developed  numerical  methods  could  very  well  solve  these  problems 
effectively.  These  numerical  methods  can  solve  not  only  range-dependent 
problems  but  also  can  handle  irregular  boundaries  with  arbitrary  boundary 
conditions.  In  this  report,  as  a start,  two  accurate  general  purpose 


‘to*  fit 


20.  (Cont'd) 

approaches  are  presented  for  the  solution  of  variable  coefficient  parabolic 
wave  equations. 

In  a finite-difference  approach,  techniques  are  derived  from  both  the 
conventional  explicit  and  implicit  schemes,  and  the  associated  convergence 
theory  is  thoroughly  analyzed.  The  techniques  are  found  to  be  general  purpose 
and  to  provide  reasonable  accuracy. 

In  an  ordinary  differential  equation  approach  the  parabolic  equation  is 
treated  as  a system  of  equations  in  which  the  second  partial  derivative  with 
respect  to  the  space  variable  is  discretized  by  means  of  a second  order  central  i 
difference  (also  known  as  the  Method  of  Lines) . Nonlinear  multistep  (NLMS)  and 
linear  multistep  (LMS)  methods  are  used  as  predictor-and-corrector  for  solving 
this  system.  A built-in  variable  step-size  technique  gives  the  desired  accuracy 
The  theory  with  regard  to  consistency,  stability,  and  convergence  has  been  veryl 
well  developed  for  both  the  NLMS  and  LMS  methods,  thus  ensuring  the  convergence f 
of  this  procedure.  | 

I 

A practical  treatment  of  irregular  boundaries  is  described  with  an 
illustrative  example,  and  a selected  set  of  experimental  numerical  results  are 
included  to  demonstrate  the  validity  of  these  approaches. 
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NUMERICAL  SOLUTIONS  OF  UNDERWATER  ACOUSTIC 
WAVE  PROPAGATION  PROBLEMS 

1.  INTRODUCTION 


Significant  contributions  to  solving  seismic  wave  propagation  problems 
have  been  made  by  applying  finite  difference  (FD)  methods however,  lit- 
tle has  been  done  in  applying  FD  methods  to  the  solution  of  underwater  acous- 
tic wave  propagation  problems.  Some  literature***^  exists,  but  further 
improvement  and  efficient  applications  of  FD  methods  to  underwater  acoustic 
wave  propagation  problems  are  required.  The  basic  problems  involved  in  the 
application  of  FD  methods  — speed,  accuracy,  and  memory  capacity  — were 
never  thoroughly  analyzed.  A study  of  the  efficiency  of  FD  methods  in 
solving  general  underwater  acoustic  wave  propagation  problems  was  never 
performed.  In  addition,  the  theory  involved  in  the  FD  methods  (which  offers 
important  information  regarding  convergence  and  error  control)  was  completely 
neglected. 

FD  methods  are  a general  purpose  scheme  and  have  very  few  restric- 
tions. ^8  a start,  we  will  search  for  efficient  solutions  of  parabolic 
wave  equations.  In  this  report  we  will: 

1.  Formulate  both  explicit  and  implicit  FD  schemes  for  variable  coef- 
ficient parabolic  equations  and  examine  whether  a more  efficient  FD 
method  can  be  developed  for  underwater  acoustic  wave  propagation 
problems, 

2.  Discuss  the  theory  of  our  FD  development  and  estimate  the  error, 

3.  Demonstrate  the  validity  of  our  FD  techniques. 

In  addition  to  the  FD  approach,  an  ordinary  differential  equation  (ODE) 
approach  was  taken  and  was  found  reasonably  effective  for  handling  parabolic 
equations.  Numerical  analysts  have  frequently  remarked  that  FD  methods  would 
be  enhanced  if  a variable  step  length  could  be  adopted  to  solve  partial  dif- 
ferential equations.  The  ODE  approach,  used  in  conjunction  with  the  FD  ap- 
proach, certainly  offers  this  advantage.  Our  method  of  attack  is  to  dis- 
cretize the  second  partial  derivative  with  respect  to  space  variables  into  FD 
representations  and  then  to  transform  the  whole  parabolic  equation  into  a 
system  of  ODEs.  The  advantage  of  using  this  approach  is  that  the  theory  of 
numerical  solutions  of  ODEs  is  very  well  developed  and  powerful  computer 
software  exists.  To  gain  speed  while  maintaining  accuracy,  we  incorporate  a 
variable  step-size  technique  with  nonlinear  multistep  (NLMS)  and  linear 
multistep  (LMS)  methods  to  handle  the  automatic  step-size  adjustment.  From 
test  results  obtained  to  date,  we  find  that  the  ODE  approach  may  turn  out 
to  be  the  most  desirable. 

In  addition  to  the  FD  material,  we  will  also 


s 
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4.  Discuss  Che  formulation  of  ODE  solutions  and  state  the  well  devel- 
oped theory, 

5.  Describe  how  to  automatically  adjust  the  step  size  by  means  of  the 
variable  step-size  technique. 

Prior  to  the  discussion  of  these  numerical  approaches,  a section  addresses 
some  of  the  theoretical  background.  After  the  descriptions  of  these 
approaches,  the  relative  merits  and  disadvantages  are  discussed.  Finally,  a 
set  of  test  results  is  given  to  demonstrate  the  validity  of  our  approaches. 
The  accurate  solution  of  parabolic  test  examples  numerically  disclose  how 
well  the  parabolic  equation  provides  a solution  of  the  sound  propagation 
problem.  A special  section  presents  a fresh,  general  treatment  of  irregular 
bottom  descriptions  with  arbitrary  boundary  conditions.  Some  conclusions  are 
provided . 

The  experimental  programs  are  written  in  ANSI  FORTRAN  language  and  have 
been  checked  out  on  the  Center's  PDP  11/70  computer.  Since  these  programs 
are  not  yet  finalized,  a listing  is  omitted.  A comprehensive  document 
describing  the  computer  model  will  be  provided  separately.  The  theoretical 
discretization  errors  are  given  in  this  report.  Since  computational  errors 
are  heavily  involved,  the  total  error  analysis  will  be  presented  in  the  com- 
puter model  documentation. 


s 
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2.  THEORETICAL  BACKGROUND 


We  begin  with  a cylindrically  synjnetric  acoustic  wave  equation  in  cylin- 
drical coordinates: 


f-f  + r 7r  + 7~2  + konZ(r»z)P  * 0 * 
dr  °z 


(2-1) 


where  p is  the  acoustic  pressure,  k0  is  the  reference  wave  number,  n(r,z) 
is  the  index  of  refraction,  and  r is  the  range  variable.  Let  p * u(r,z)v(r) 
where  v(r)  is  strongly  dependent  on  r,  but  u(r,z)  is  only  weakly  dependent 


on  r. 


By  substituting  p ■ uv  into  equation  (2-1),  rearranging  the  terms,  and 
2 

using  k as  a separation  constant,  we  obtain  the  following  reflected  field 
and  transmitted  field  parabolic  wave  equations: 


1 2 
v + — v + k v 
rr  r r o 


(2-2) 


u + 
rr 


* (l  + 2 3v\ 
zz  \r  V r 


+ k^  (n^-l)u  * 0. 
o 


(2-3) 


After  we  introduce  the  far-field  approximation  (k0r  y>  1),  the  exact  solu- 
tion of  equation  (2-2)  is  known: 


v(r)  - H(1)(k 
o o 


r>*^k*ltk#r'5>  • 


(2-4) 


Using  equation  (2-4)  to  eliminate  v from  equation  (2-3),  we  obtain 


u + u + 2ik  u + k^(n^'-l)u  ■ 0. 
rr  zz  o r o 


(2-5) 


If  we  assume  that  inhomogeneities  vary  slowly  with  range,  the  reflected 
field  can  be  neglected,  which  leads  to  the  following  parabolic  wave  equation: 


ik  (n  -1) 

o 1 

u ■ , u ♦ t: — u . 

r 2 2k  zz 

o 


(2-6) 


Performing  parabolic  approximations  to  the  reduced  wave  equation  will  result 
in  various  slightly  different  forms . 2, 5 , 6 , 10  choose  to  deal  with  equa- 
tion (2-6)  because  it  is  commonly  recognized  in  the  underwater  acoustics 
community,  and  existing  PE  (parabolic  equation)  models,  such  as  the  split- 
stepH,  are  directed  toward  the  solution  of  equation  (2-6). 
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Note  that  in  the  shallow  water  environment  or  wherever  bottom  interac- 
tion is  important,  the  exact  boundary  condition  must  be  satisfied;  therefore, 
more  general  purpose  methods  are  sought  for  the  solution  of  equation  (2-6). 
Existing  methods  using  the  fast  Fourier  transformation  (FFT)  are  inappropri- 
ate for  such  problems  because  these  methods  reflect  the  entire  medium  across 
the  surface;  however,  below  the  physical  ocean,  an  absorbing  layer  is  intro- 
duced, permitting  the  infinite  transforms  to  be  truncated,  so  that  the  FFT 
becomes  applicable.  The  numerical  ODE  and  the  FD  approaches  are  advantageous 
because  they  do  not  introduce  an  artificial  bottom. 
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3.  NUMERICAL  SOLUTIONS 

3.1  FINITE  DIFFERENCE  APPROACH 

3.1.1  Formulation 


We  consider 

u ■ a(k  ,r,z)u  + b(k  ,r,z)u  ■ Lu.  (3-1) 

r o o’  zz 

Let  D be  the  FD  operator,  and 

5Z  be  the  central-difference  operator  in  the  z-direction. 

D and  6Z  are  related  by 


D 


2 . 

— sin 


-1 


Where  h * Az,  we  use  k * Ar.  Therefore, 


D 


2 


1*  52  ♦ ^ 64 

12  z 90  z 


(3-2) 

!i 


We  choose  to  write 


D 


2 


[1 


p(k,h)] 


(3-3) 


p(k,h)  is  a parameter  to  be  determined  such  that  the  choice  of  p(k,h)  will 
minimize  the  initial  local  discretization  error  for  implicit  schemes. 

Further,  we  assume  that  p(k,h)  + 1.  If  p(k,h)  * 1,  equation  (3-1) 
reduces  to  ur  * au,  which  is  trivial.  In  fact,  we  will  demonstrate  that 


p(k,h)  ■ 0 gives  the  Crank-Nicolson  formula, 
p(k,h)  - h2/6k  gives  the  Douglas  formula. 

We  will  make  an  attempt  to  find  a p(k,h)  such  that  the  order  of  the  error  is 
the  smallest,  if  possible. 


3. 1.1.1  Explicit  Methods.  We  start  with  the  formulation  of  explicit 
methods.  Using  the  Taylor  expansion,  we  can  obtain  a two-level  scheme  for 
equation  (3-1): 


u(r  + Ar,z) 


k 


1 + k^—  + — k^ 

1 V 2!  K 


■i,* 

3r 


-■) 


u(r,z) 


k -a. 

e 3r  u(r,z) 


(3-4) 
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If  we  let  z * mh,  r 3 nk,  and  u(r,z)  ■ u(nk,mh)  * u”,  then  equation  (3-4)  can 
represent  an  explicit  scheme: 


n+1 

u 

m 


(3-5) 


To  solve  equation  (3-1),  we  could  use  equation  (3-5)  and  retain  only  the 
second  order  difference;  the  explicit  formula  becomes 


+ k h 


r ) Um  * (X  + ak  k 5z)Um  * 


(3-6) 


In  general,  we  drop  the  superscripts  and  subscripts  of  a and  b for  econ- 
omy in  writing;  however,  in  the  critical  places,  we  will  use  the  indexes  to 
make  our  formulas  clearly  readable.  Using  the  second  order  central  differ- 
2 

ence  for  5 , we  find  that 
z' 


(1  + ak) 


n b / n n n \ 

m 2 ym+1  m m-1  / 


(3-7) 


If  we  write  1 + ak  * otn,  k ■ 0n,  we  can  express  equation  (3-7)  in  a matrix 

TO  , Z TO 

n 

form: 


ods,  we  use 
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To  solve  equation  (3-1),  we  u*e 


[l  - -k(a  + bD2)]  un+1  * [l  + -k(a  + bD2)  u11. 
l m z m 


(3-10) 


Again,  a second  order  central  difference  is  used  for  D2  in  equation  (3-10); 
we  obtain 

fs2  1 r 62  1 

1 - ka  - |kb  -f  (1-P(k,h))  un+1  «'l  + ka  ♦ kb  -f  (1-P(k,h))  u"  . 

2 2 Lh2  J m 2 2 Lh2  J m 


Simplifying  the  above,  we  find  that 


= 1 - -jka  + bs(l-P)  u"+1  - jbs(l-P)  (n^\  + u£jJ 

m 1 + ka  - bs(l-P)  un  + ks(l-P)  /u  + u11 

z m z v m+ 1 m- 1 1 * 


where  s - ^ • 
If  we  write 


a = 1 - ka,  B = bs(l-P),  P = P(k,h),  Y = 1 + yka  , 


(3-11) 


(3-12) 


(3-13) 


x - xn+1  = ctn+1  + Bn+1  - 1 - kan+1  + bn+1  s(l-P), 
mm  m m Z m m 


Y * Yn  = Yn  - Bn  =»  1 + ka11  - bn  s(l-P), 

m m m m Z m m 

and  equate  both  sides  (LHS  ■ RHS),  we  find  that 


1 on+1  n+1  . 
- — P u , + 

2 m m+1 


(^n+1  ftn+l\  n+1  1 fin+l  n+1 

a + P \u  - -r  P u , 
m m J m 2 m m-1 


. 1 0n  un  , + kn  - Qn\  u"  + jBnun. 
2 m m+1  l m m J m 2 m m-1 

Equation  (3-14)  can  be  written  in  matrix  form  as: 


(3-14) 
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The  two  components  of  the  first  column  vector  on  the  RKS  are  two  boundary 
points,  and  the  two  components  of  the  last  column  vector  on  the  RHS  are  two 
boundary  points  on  the  initial  line.  (Note;  If  we  select  P * 0,  equation 
(3-15)  reduces  to  the  Crank-Nicolson  scheme;  if  we  select  P * h^/6k,  equa- 
tion (3-15)  reduces  to  the  Douglas  scheme.)  At  this  point,  we  cannot  yet 
make  a choice  of  P to  arrive  at  a new  method.  We  defer  this  until  after  we 
have  developed  the  consistency  of  the  method. 


3.1.2  Consistency.  The  conventional  definition  of  consistency  states 
that  an  FD  approximation  to  a parabolic  equation  is  consistent  if 


Truncation  error 
k 


0 as  h,  k 0. 


—r 
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However,  we  adopt  the  definition  of  consistency  in  the  sense  of  Keller^  as 
follows : 

Let 

£°[u;h,k]  be  an  FD  approximation  to  equation  (3-1), 

£ [u]  be  the  true  operator. 

Define 

T [u;h,k]  - £[U]  - £°[u;h,k] . 

If 


lim  T [u;h,k]  0, 

h,k  0 

we  say  that  the  method  is  consistent,  meaning  that  the  FD  operator  is  con- 
sistent with  the  true  operator. 

How  we  proceed  to  develop  the  concept  of  consistency  and  to  obtain  the 
"initial  local  discretization  error"  (usually  called  the  "truncation  error"). 

Expanding  un+^,  un  , , un  , upon  u11,  using  the  Taylor  expansion,  and  sub- 
m m+ 1 m- 1 m 

stituting  the  results  into  equation  (3-7);  we  obtain 

n+1  /iJ.,\n  b./n  ,,  n . n \ 

u - (l+ak)u r-k  (u  . - 2u  + u . 

m m 2 \ m+1  m m-1  J 


m ' ' m ' tti 


• • • 


(3-16) 


The  terms  inside  the  { } of  equation  (3-16)  » 0 because  they  satisfy 

equation  (3-1).  Let  E[e]  indicate  the  principal  part  of  the  initial  local 
discretization  error.  Then,  we  have 

E[e]  * 0(k^  + kh^) . 

As  hjk-*’  0,  limT[u;h,k]  0;  therefore,  method  (3-7)  is  consistent. 
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..  n+1  n+1  n n n . ..  _ , 

Similarly,  expanding  u_., , u_  .,  u. , , u_  . upon  u , using  the  Taylor 

m+  i m- 1 m+  l m—  1 m 

expansion,  substituting  the  results  into  equation  (3-14),  and  simplifying,  we 
obtain 

j 1 Qn+1.  n+1  . /n+1  . Qn+1\  n+1  1 an+l  n+1 

l 2 m m+1  \m  m/m  2m  m-1 


1 1 Qn  n 
2 m um+l 


n Qn\  n^lQn  n 
W - S ) u +tp  u . 
Vm  m/m  2 m m-1 


n /Su'' 
-kau  + k 
m 


(&) 

m 


2 l 3ju\  ” 


- bsh  , 
V3zZ 


m 


. [-j-®;- •*(£»);■ -£(&)] 

. . 2 /32u\  n h2.  / 33u  \ n h2  k2/  3 4u  \ n 

* pb,b  (i?)  * pb5^  fe)  „ * pbs^  irfej 

m m 


b,a-4  ($) " • 

tn 


• • • • 


(3-17) 


The  terms  inside  the  { } of  the  RHS  of  equation  (3-17)  * 0 because  they 

satisfy  equation  (3-1).  If  a and  b are  range  independent,  we  see  that  the 
terms  inside  the  £ 3 and  the  ( ) of  the  RHS  of  equation  (3-17)  all  * 0 

because  the  terms  inside  the  [ } satisfy  urr  - aur  - buzzr  * 0 and  the 
terms  inside  the  ( ) satisfy  u^j.  - au^  - buzzrr  * 0.  Let  E[I]  indi- 

cate the  principal  part  of  the  initial  local  discretization  error  of  the  im- 
plicit scheme  of  equation  (3-14).  ^ — 


For  p*0, 


E [ I ] ■ 0(k2  + kh2),  a Crank-Nicolson  error; 


p • 0(kh2),  E[I]  * 0(k3  + kh^),  a Douglas  error. 

It  does  not  seem  likely  that  a p can  be  chosen  such  that  E[I]  is  smaller  than 
the  Douglas  error.  However,  we  see  that  lim  x [u;h,k]  ■+  0 ; therefore,  the 

h,k+0 

implicit  methods  of  formula  (3-15)  are  consistent  for  range- independent  coef- 
ficients a and  b.  For  the  range-dependent  case,  the  terns  inside  the  [ 3 
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and  Che  ( ) of  Che  RHS  of  equation  (3-17)  do  noc  vanish,  so  ChaC  we  will  ob- 

cain  a larger  E[I].  For  all  p's,  we  can  see  ChaC 

E [I]  - 6(k2  ♦ kh*  ) . 

s 

However,  lim  x[u;h,k]  ■*  0;  Cherefore,  Che  implicic  mechods  of  formula  (3-15) 
h,VwO 

are  consisCenC  also  for  range-dependenc  coef f iciencs . 

3.1.3  Stability.  The  general  concepc  of  scabilicy  scaces  chat  the 
difference  becween  Che  CheoreCical  and  numerical  solucions  remains  bounded  as 
Che  range  sCep  n increases,  provided  Che  range  increment  k remains  fixed  for 
all  space  seeps  m.  To  find  ouC  whether  a method  is  stable  or  noc,  we  examine 
Che  satisfaction  of  Che  scabilicy  condition;  this  condition  can  be  derived  by 
means  of  familiar  methods  such  as  Von  Neumann's,  Che  macrix,  or  Che  Fourier 
series . 

We  shall  first  derive  Che  scabilicy  condition  for  Che  explicit  mechods, 
equation  (3-7).  We  apply  Von  Neixsann's  criterion  of  stability  to  equation 
(3-7)  by  seeking  a solution  in  the  form  ear  e1^1.  Substituting  the  solu- 
tion into  equation  (3-7),  we  get 

_n+l  iojmh  _n/t  n,  n . iujh 

S e - F (1  ♦ a*k  - 2b  s)e 

ED  ID 


where  £ * e2^. 

Simplification  of  equation  (3-18)  gives 

_ , A n,  . . n . 2 /wh\ 

E • 1 + a k - 4s  b sin  {-= — I . 
s m m \2  / 

|^|  < 1 is  required  to  give  the  stability  condition;  i.e., 

I 1 + - 4s  b”  sin2  (f)\  < 1.  (3-19) 

When  we  use  formula  (3-7)  to  solve  the  example  equation,  ur  ■ u22,  we 
arrive  at  Mitchell's^  results; 


♦ 5n  bn 

» in 


i 

e 


iu)(m+l)h  ^ iui(m-l)h 


♦ e 


— 

i 

» 


(3-18) 


‘ 1 £ 1 * 4\  sin2  (a?)  i 1 » 

h 

. . . k 1 

which  implies  s * —?  <_  — for  stability. 

h 

In  general,  a(k0,r,z)  and  b(kQ,r,z)  are  functions  of  r and  z;  this 
requires  a thorough  examination  of  condition  (3-19).  In  the  case 
a(k0,r,z)  • 0,  we  need 
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As  long  as  s,  b° 
isfied.  The  LH? 


> o, 


1 < 1 - 4sb“  sin2  (Sjfe)  < ♦ 1 . (3-20) 

the  RHS  inequality  of  equation  (3-20)  is  trivially  sat- 


gives 


In  the  case  where  bn 

m 


is  purely  imaginary,  we  need 


| 1 - 4sb”  sin2  (I*)  | < 1.  (3-22) 

Condition  (3-22)  does  not  hold  for  s > 0,  bn  purely  imaginary;  therefore, 

in 

formula  (3-7)  is  not  stable  for  problems  with  zero  coefficient  a(k0,r,z) 
and  a purely  imaginary  b(kQ,r,z). 

When  a(kQ,r,z)  ^ 0,  we  need  inequality  (3-19)  to  hold.  In  our  appli- 
cations, a(kQ,r,z)  and  b(k0,r,z)  are  both  purely  imaginary.  Let 
a » iaR,  b * ib^;  we  then  have 

|l  + i faRk  - 4sbR  sin2  | ,<  1.  (3-23) 

L 

Obviously,  inequality  (3-23)  does  not  hold.  To  hold  the  inequality,  we  must 
have 

||  V - 4,bR sir,J  (r) ' °> 

which  implies  that  the  condition  of  stability  is 

h2  - 4 ^ sin2  (yl)  . (3-24) 

Equation  (3-24)  holds  for  h • 0,  which  implies  the  instability  of  scheme 
(3-7).  To  overcome  this  difficulty,  a parameter  X can  be  introduced  to 
obtain  an  equivalent  equation  of  scheme  (3-7);  namely, 

vr  " [*<Vr’z)  * Xlv  + b(ko,r,z)vzz’  (3-7)' 

where  v ■ e”*r  u.  Then,  a similar  inequality  of  equation  (3-19)  is 
obtained;  namely, 
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| (1+X)  ♦ i ^a^k  - 4sbRsin2  I < 


(3-19)' 


The  values  k,  h,  and  X can  be  chosen  such  that  equation  (3-19)'  is  satisfied; 
therefore,  scheme  (3-7)'  is  conditionally  stable. 

Now  let  us  turn  to  the  stability  of  implicit  methods,  equation  (3-14). 
When  a(k0,r,z)  ■ 0,  b(k0,r,z)  • 1,  and  0 ■ 0;  equation  (3-14)  reduces  to 
the  example  equation  Uy  * uzz  and  our  formula  reduces  to  the  Crank- 
Nicolson  formula.  Using  Von  Neunann's  method  again,  we  find  that  the  stabil- 
ity condition  is 


, „ 1 - s .1  - cos(a)h)]  . . 

1 - ITT  1 - cosUfi  i l* 


(3-25) 


We  see  that  condition  (3-25)  is  satisfied  for  a > 0.  In  fact,  s ■ k/h2  is 
always  > 0;  therefore,  we  see  that  the  Crank-Nicolson  formula  is  uncondition- 
ally stable  for  solvipg  ur  » uzz. 

However,  for  the  equation  with  complex  coefficients,  we  cannot  take  the 
unconditional  stability  for  granted  when  we  are  using  the  generalized  Crank- 
Nicolson  formula;  the  stability  of  formula  (3-14)  needs  a thorough  investiga- 
tion. We  shall  apply  Von  Neumann's  method  to  formula  (3-14)  to  derive  a con- 
dition and  examine  the  condition  in  detail. 

Define  £ * eak  e^*1.  Substituting  C into  formula  (3-14),  we  get 

b a(n+l)k  iw(m+l)h  . _ ct(n+l)k  iwmh  b a(n+l)k  iw(m-l)h 
- - se  e + Ae  e • se  e 

b ank  iu>(m+l)h  ^ „ ank  iwmh  ^ b ank  iw(m-l)h 

• -r  se  e + Ye  e ♦ -r  se  e 


Simplifying,  we  find  that 


Y ♦ ^ s [2  cos(wh)] 
X - ^ s [2  cos(wh)3 


(3-26) 


which  gives  the  condition 


1 - bs  [ 

1 - cos(wh)]  - ak 

1 ♦ bs  [ 

1 - cos(ojh)]  ♦ ak 

< l . 


(3-27) 


Where  a(k0,r,z),  b(k0r,z)  > 0,  inequality  (3-27)  is  satisfied  for  all 
s > 0;  therefore,  the  unconditional  stability  obtains.  In  our  applications, 
a(k0,r,z)  and  b(ko,r,z)  are  both  purely  imaginary.  Specifically,  we  have 
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1 - ~ s [l  - cos  (wh)J  - J 

0 

! 2 

n2(r,z)  - l| 

1 + s Fl  — cos(wh)J  + ^ 

0 

ik 

0 

2 

n2(r,z)  - lj 

1 1 • 


(3-28) 


Equation  (3-28)  can  be  expressed  in  the  short  form 

1 - iX 


1 + iX 


< 1, 


(3-29) 


where 


X - [l  - cos(u)h)]  + [n2(r,z)  - l]  . 

o 

From  equation  (3-29)  we  see  that 

always  * 1 for  all  real  X . 

Therefore,  the  FD  scheme,  equation  (3-14),  is  stable  when  it  is  applied  to 
solve  our  parabolic  wave  equation. 

3.1.4  Convergence.  We  have  developed  the  consistency  and  the  stabil- 
ity of  both  explicit  and  implicit  FD  schemes,  formulas  (3-7)  and  (3—14).  Now 
we  want  to  show  that  our  FD  schemes  are  convergent.  Let 

"t.s."  stand  for  the  theoretical  solution  of  equation  (3-1), 

"n.s."  stand  for  the  numerical  solution  of  equation  (3-1), 

"f.s."  stand  for  the  finite-difference  solution  of  equation  (3-1). 


1 - iX 
1 + iX 


The  norm  inequality  shows  that 

1 1 t.s.  - n.s. ||  < 1 1 t.s.  - f.s. ||  + 1 1 f.s.  - n.s. ||  . 

Note  that  by  applying  the  consistency  to  the  fi.rst  norm  of  the  RHS  and  by 
applying  the  stability  with  the  error  control  to  the  second  norm  of  the  RHS, 
we  have  established  the  convergence. 

3.1.5  Discretization  Errors.  As  a result  of  their  consistency,  the 
local  initial  discretization  errors  of  explicit  and  implicit  FD  methods  are 
given  by  formulas  (3-16)  and  (3-17),  whose  principal  parts  possess  the  fol- 
lowing expressions: 
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Discretize  the  uzz  portion  by  a second  order  central  difference: 


15 


Decompose  Che  above  matrix  form  inCo 


which  is  in  Che  form 

u'  ■ Au  + g(r,z,u)  (3-38) 

3.2.2  Methods  Under  Consideration  (NLMS  and  LMS) 

The  above  form  (equation  (3-38))  is  the  equation  to  which  the  NLMS  meth- 
ods are  applicable.  To  simplify  Che  application,  we  apply  an  NLM-l-step  pre- 
dictor and  a corrector  with  a built-in  variable  step  size. 


As  a predictor,  the  solution  is 
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n+1  _ Ah  n . . r Aho  n,  n. 

i * e u + h lj(Ah)  (I-e  )J  g (r  ,z,u  ). 


(3-39) 


As  a corrector,  the  solution  is 


n+1  Ah 
u ■ e 


un  + h jj(Ah)  ^ 


I + (I-Ah)eAh  gn 


+ (I  + Ah  - eAh)  gn+^(rn+^,z,un+^)  . (3-40) 


A PC3  procedure  is  built-in;  i.e.,  the  procedure  predicts  and  corrects  at 
most  three  times. 


NLMS  methods  are  designed  to  be  effective  in  solving  equations  whose  g 
is  either  slowly  varying  in  r or  is  a low  order  polynomial  in  r and  whose 
eigenvalues  of  A have  negative  real  parts  and  differ  greatly  in  magnitude. 
The  selection  of  step-size  h can  be  made  approximately  by 


V<Ah> 


h < 1, 


where  has  been  defined  in  formula  (3-40). 


When  g is  a constant,  any  h can  satisfy  the  above  inequality;  therefore, 
NLMS  methods  allow  the  use  of  a large  step  size.  In  this  case,  an  accurate 
computation  of  e^*  will  give  accurate  results  with  fast  speed.  Therefore, 
NLMS  methods  have  an  excellent  application  to  underwater  acoustic  .range- 
independent  problems  with  plane  parallel  boundary  conditions  where  the  field 
vanishes  on  both  boundaries. 

3.2.3  Stability,  Consistency,  and  Convergence 


The  theory  with  respect  to  consistency,  stability,  and  convergence  has 
been  well  developed. 13  We  summarize  the  theory  below. 


NLMS  methods  take  the  expression 


k 

Z 

i-0 


a. 

i 


Ah(k-i) 
e u 


n+i 


h 


k 

Z 

i-0 


<t>ki(Ah)8 


n+i  * 


(3-41) 


The  characteristic  polynomial  of  NLMS  methods  is 

p(X,C)  - eXAh  Z ct.C1  . (3-42) 

i-0  1 

Since  we  select  k-1,  ct^-l,  and  ct^-i  - -1,  the  root  of  p(X,C)  is  1 and 
simple;  therefore,  the  NLMS  methods  are  stable. 

The  consistency  condition  is  self-contained  because  NLMS  methods  are 
formulated  to  yield  consistency.  Then,  by  the  convergence  theorem,  stabil- 
ity + consistency  *+  convergence. 
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We  also  tried  LMS  methods,  which  can  be  simply  obtained  from  NLMS  meth- 
ods by  letting  | | A|  | -*•  0 so  that  the  theory  is  automatically  applicable. 

The  predictor  we  used  is 

un+l  ■ un  + hfn(r,z,un),  (3-43) 

which  is  known  as  the  Adams-Bashf orth  method  and  also  carries  the  name  Euler 
method. 

The  corrector  we  used  is 

n+1  n h r,n  ,n+l,  n+1.*]  ,, 

u » u ♦ — If  + f (r,z,u  )J  , (3-44) 

which  is  known  as  the  Adams-Moul ton  method  or  the  trapezoidal  rule,  f here 
is  the  RHS  of  equation  (3-36). 

3.2.4  Discretization  Errors 

The  initial  local  discretization  errors  of  LMS  and  NLMS  methods  have 
been  worked  out  in  detail^»14.  The  error  terms  are  listed  below  in 
relation  to  the  methods  we  have  applied. 

E [ AB]  ■ Adams-Bashf orth  error 

- hp+2  u(p+2)  (?)  [(-1)P+1  J*  dsj  - 0(hP+2)  (3-45) 

E[AM]  * Adams-Moul ton  error 

(£)  [(-l)'*'  J.;  (j:,)  J>]  - 0(h**h  (3-46) 

E [NLMS]  ■ Nonlinear  multistep  error 

» (constant)  | | g^P  ^(r,u)||hP  2 * 0(hP  2),  (3-47) 

where  p " order  of  the  method, 


3.2.5  Variable  Step-Size  Technique  and  Error  Controls 

Definitions:  PCm  stands  for  "predict-and-correct-m-times ," 

e is  the  user  required  tolerance. 
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A variable  step-size  employs  the  PCm  procedure  to  satisfy  the  user's 
required  tolerance. 


Applying  a predictor  (explicit  method),  one  can  use 

k-1 


n+k 


i r k_i 

j-  - z «■ 

a k i*0 


^(k-i)^  h Z ♦ (Ah)g 

i*0 


where  u^+k  is  expressed  in  terms  of  previous  un+£  values. 

If  a corrector  (implicit  method)  is  used,  one  must  have 


n+k 

which  is  of  the  form 

where 


u * G(u) , 


u * «n+k* 

The  successive  iterative  form  gives 

U(V+1)  „ G(u(v)) 

for  any  initial  vector  u^O) . 


(3-48) 


[-  Z a.eAh(k_l)u  +h  I <t>  .(Ah)g  1 , 

ak  L i«0  1 0+1  i*0  kl  n+1J 


(3-49) 


(3-50) 


(3-51) 


Let  G(u)  be  defined  for  1 1 u)  | < ®,  and  let  there  exist  a constant  K 
such  that  0 < K < 1.  Then  G(u)  satisfies  the  condition 


| | G(u*)  - G(u)|  | < K | | u*  - u| 


(3-52) 


Using  the  definition  of  G(u)  — formula  (3-49)  — and  the  fact  that  g(r,z,u) 
satisfies  the  Lipschitz  condition  with  a Lipschitz  constant  L,  one  can  see 
that  equation  (3-52)  is  satisfied  by 


h||*k  k(Ah)M 

K - ^ L . 


(3-53) 


For  the  iterative  procedure,  equation  (3-51),  to  converge  for  arbitrary 
u^O) , K is  required  to  be  less  than  Is 

h||4>kk(Ah)|| 


K < 1 + 


L < 1 . 


(3-54) 
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This  is  the  condition  of  the  corrector's  convergence.  Conventionally,  we 
select  oik  ■ 1 for  computational  convenience. 

Now  we  describe  how  we  develop  the  variable  step-size  technique.  If  a 
predictor  is  accurate,  only  one  correction  is  needed.  In  our  implementation 
we  limit  m£  3;  thus,  we  have  a PC^  procedure,  which  can  be  described  by 
the  following  diagram. 


© , the  latest  un+^,  meets  the  corrector's  convergence;  therefore,  un+^ 

is  the  solution.  When  the  corrector's  convergence  is  met,  we  double  the  step 
size  and  go  to  (§)  to  continue  the  solution.  When  the  third  try  fails  to 
meet  the  convergence,  we  halve  the  step  size  and  go  to  (§)  to  restart  the 
solution.  Since  there  exists  an  h that  leads  to  the  corrector's  convergence, 
this  procedure  is  a convergent  procedure;  however,  when  h is  extremely  small 
such  that  t+h  * h in  the  machine,  the  program  will  provide  a message  and  halt 
the  computation. 

When  it  is  less  than  the  user  required  tolerance,  the  accuracy  is  con- 
trolled by  examining  the  norm  between  two  consecutive  computations. 
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4.  A NUMERICAL  TREATMENT  OF  BOUNDARIES 

Consider  the  solution  of  the  parabolic  wave  equation  of  the  form 

u = a(k  ,r,z)u  + b(k  ,r,z)u  . (4-1) 

r o o zz 

Assume  that  the  numerical  methods  (FD  and  ODE)  are  to  be  used  to  solve  the 
problem  in  the  rectangular  region: 


BOUNDARY 
P01NT-*o — 


SURFACE  BOUNDARY 


INITIAL - 
VALUES 


• * 


BOTTOM  BOUNDARY 


BOUNDARY 

POINT 


Our  methods  require  the  initial  conditions  and  two  boundary  points 
before  we  can  proceed  to  find  the  solution  at  the  next  range  line.  With 
this  input  information,  we  can  classify  the  boundary  descriptions  for  three 
different  cases,  and  discuss  the  treatment  of  each  case  separately. 

We  further  assume  that  the  boundary  condition  is  a Neumann  condition 
such  that  p^  = 0,  where  p is  the  solution  of  V^p  + k^p  = 0, 


CASE  1 : FLAT  BOTTOM  BOUNDARY 


In  this  case, 


BOTTOM  BOUNDARY 
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Since  p * H^^(kr)u,  p * 0 *►  u =0. 
r o z z 

n+ 1 n+ 1 

u . - u 

To  find  un+,  , we  see  that  t — — — = 0;  therefore, 

m+1  tiz 


CASE  2:  SLOPING  BOTTOM  (SHALLOW  TO  DEEP  WATER) 
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We  find  un+}  as  follows: 
m+1 


^ cosa  - sinOt  » 0 at  u 


n+1 

^m+l 


(4-4) 


n+1 


Also  at  u s u , , u must  satisfy  parabolic  equation  (4-1);  therefore, 
m+I 

equation  (4-4)  becomes 

H^(kr)u  cosa  - (H^(kr))  u + H^(kr)u  sina  * 0 . 

o z [_  o r o rj 


And  therefore, 


^ (kr)cosdu 
o z 


- [Hl°(kr)]r 


sindu  - ^ (kr) sind( au+bu  ) = 0. 

o zz 


Simplifying,  we  obtain 


cota 

u - — r—  u + 
zz  b z 


(H^n(kr))r 


H(1)(kr) 
L o 


+ a 


b u " °» 


(4-5) 


which  is  a second-order  scalar  ODE. 


Knowing  6z,  un+\  we  can  find 
in 


n+1  n+1 


/ 3u\n+1 

\ 


from 


n~r  i iv  i 

u - u . /*  \ n+1 

m m-1  I du  l 


(4-6) 


Then  we  can  use  numerical  methods  to  solve  equation  (4-5).  We  make  the  fol- 
lowing transformation,  using  the  ODE  package  already  incorporated  in  the  sys- 
tem. Let 

cota 

P1  " ’ 


(R(1)(kr)) 
o r . 

— m 3 

HU;(kr) 

o 


rb  . 
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If  Pi(z)  and  P2(z)  are  constants  in  z,  we  can  solve  equation  (4-5) 
exactly;  i.e., 


The  solution  is 


u-  V(f*2  yp?-‘p2)'  ♦ c2e(r'2  M - 


(4-7) 


If  we  use  u(z0)  » u0,  u (z0)  3 u , we  can  determine  cj, 

n+l 

C2  exactly.  Then  the  solution  u can  be  calculated  exactly,  and  u . can  be 
solved.  m 1 


Once  u;, | is  found,  we  can  proceed  if  we  use  proper  care.  Note  that 

if  the  distance  between  u^+*  and  u^j  is  not  * 2(5z),  we  cannot  proceed  with 

the  present  method  because  the  depth  partition  at  this  time  level  is  not  uni- 
form. A special  treatment  is  needed  at  some  specific  advanced  time  levels. 

To  handle  this,  we  label  u"*}  as  and  define  u”T}  as  the  mid-point 

between  u and  u We  determine  un  , in  such  a way  that  the 

m m+2  m+1 
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distance  between  u^  ^ and  u^+2  is  2(6z).  This  treatment  ensures  a uniform 

partition  in  the  depth  direction.  The  u^j  can  be  found  the  same  way  by 
solving  equation  (4-5). 

A question  arises:  How  far  can  we  proceed  in  determining  un^|?  To  ease 

m+ 1 

the  programming  requirements  and  to  maintain  accuracy,  we  determine  Ar  by 


tga  •*  Ar  ■ 5z  • cota. 


The  treatment  of  case  2 can  be  applied  to  case  3.  The  only  differences 
are  a and  5z. 
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5.  EXPERIMENTAL  NUMERICAL  RESULTS 


To  implement  these  numerical  methods,  experimental  computer  programs 
have  been  developed  in  ANSI  FORTRAN  language.  The  computations  were  carried 
out  on  the  Acoustic  and  Environmental  Research  Division's  (Code  312)  PDP 
11/70  computer.  Since  these  programs  are  not  yet  finalized,  we  will  not 
describe  them  in  this  report.  A separate  document  will  be  issued  to  describe 
the  computer  model  in  detail. 

A number  of  parabolic  problems  whose  solutions  are  known  were  used  to 
test  our  numerical  methods;  the  numerical  results  show  a reasonable  accuracy. 
Three  of  the  test  problems  are  included  in  this  report  to  demonstrate  the  va- 
lidity of  the  methods.  Note  that  problem  1 (used  to  test  accuracy)  does 
not  have  any  physical  significance. 

We  generally  classify  underwater  acoustic  wave  propagation  problems  in 
two  categories:  range-independent  problems  and  range-dependent  problems. 
Range- independent  problems  have  coefficients  and  plane  parallel  boundary 
conditions  independent  of  range.  Range-dependent  problems  include  either 
irregular  bottom  boundaries  or  plane  parallel  boundaries  with  coefficients 
and/or  boundary  conditions  dependent  on  range. 

PROBLEM  1:  ACCURACY  TEST 
Equation:  Uj.  * (r^-z)u  + uzz  . 

Region  of  consideration: 


0<  r<  1 
0<Z  <7T 


Conditions:  u0  ■ u(o,z)  ■ 1 
u(r,o)  ■ 1 


u(r,»)  ■ e"Wr 
Exact  solution:  u • e"rz. 
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Numerical  Solutions: 

Results  are  tabulated  at  a range  of  1 mile. 


Depth 

(mile) 

ODE 

Implicit 

FD 

Exact 

Solution 

mm 

0.81884402 

0.81880862 

0.81873077 

mm 

0.60672331 

0.60663551 

0.60653067 

0.8 

0.44953302 

0.44939151 

0.44932896 

time 

1®229 

lm178 

A variable  step-size  technique  is  used  for  the  ODE  method;  the  user  re- 
quired tolerance  is  of  the  order  h"^.  The  FD  method  does  not  have  the 
variable  step-size  capability;  it  uses  the  optimal  step-size  (h*0.001)  deter- 
mined by  the  ODE  program. 


PROBLEM  2:  A SHALLOW  WATER  WAVE  PROPAGATION  IN  A RECTANGULAR  REGION  WITH  A 
RIGID  BOTTOM15 


Equation: 


ik 


u ■ — r—  (n  -l)u  + -z r—  u ; k • . 

r 2 2k  zz’  o c 

o o 


Region  of  consideration: 


BOTTOM 

Input  parameters  and  initial  boundary  conditions: 

Initial  field  values:  supplied  by  shallow  water  model15 
Surface  condition:  u(r,0)  ■ 0 
Bottom  boundary  condition:  ug  ■ 0 
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Source  depth:  32  ft 
Initial  sound  speed  * c0  ■ 4950  ft/ sec 
Bottom  depth  ■ 64  ft 
Frequency  * 25  Hz 
Sound  speed  profile: 


z (ft) 

SSP  (ft/sec) 

0 

4950 

64 

5000 

>64 

90000 

Nimerical  solutions:  Problem  started  at  range  * 32  ft  and  terminated  at 

range  * 1313  ft.  Three  solutions  were  obtained;  a 
normal  mode  15,  a variable  step-- size  ODE,  and  an 
explicit  FD.  The  FD  program  used  a step-size  10  times 
smaller  than  that  used  for  ODE.  The  FD  solution  was 
obtained  using  double- precision  arithmetic. 

Graphical  outputs:  Three  numerical  solutions  are  plotted  on  the  graph 

below:  depth  (ft)  versus  propagation  loss  (dB). 
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PROBLEM  3s  SLOPING  BOTTOM  (SHALLOW  TO  DEEP  WATER  ENVIRONMENT) 16 


Problem 

Background: 


Equation: 


The  solution  for  the  acoustic  field  in  a homogeneous  medium 
bounded  by  a wedge  can  be  obtained  by  the  method  of  images . 

The  solution  in  mathematical  form  and  a computer  program 
(WEDGE)  to  evaluate  the  exact  solution  are  described  in 
reference  16.  A homogeneous  medium,  characterized  by  a sound  speed 
c0,  bounded  by  a horizontal  surface  and  sloping  bottom  that 
makes  an  angle  of  45  degrees  with  the  horizontal  surface  is 
described  below  under  Region  of  consideration. 

The  parabolic  equation  represents  the  wave  propagation  in  the 
range  direction  after  the  parabolic  decomposition  is  found  to  be 


2kg Uzz 


Region  of  consideration: 


RANGE 


RECEIVER 
LINE 


A point  source  is  placed  50  ft  below  the  surface,  and  the  receiver  is  located 
100  ft  from  the  source.  The  depth  of  the  wedge  at  the  source  is  100  ft,  the 
frequency  equals  80  Hz,  and  c0  equals  5000  ft/sec.  The  exact  solution 
above  indicates  the  solution  of  the  acoustic  wave  equation.  We  attempt  to 
find  the  parabolic  equation  solution  u in  the  shaded  region,  and  then  we 

attempt  to  compare  H^^(kr)u  against  the  acoustic  wave  equation 

solution  computed  by  the  WEDGE  program. 
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Conditions:  Initial  u(r0,x)  values  are  supplied  by  the  WEDGE  program. 
u(r,o)  • 0 

u(r, bottom)  ■ ujj  ■ 0. 

Exact  solutions: 

The  exact  solution  of  the  acoustic  wave  equation  is  obtained  from  the 
WEDGE  program. 

The  exact  solution  of  the  parabolic  equation  is  not  known. 

Numerical  solutions: 

Numerical  solutions  are  produced  by  the  variable  step-sire  ODE  method. 
The  Neumann  bottom  condition  was  treated  by  formula  (4-5).  Solutions  of  the 

parabolic  equation  are  multiplied  by  H^^(kr)  to  give  an  approximate 

solution  of  the  acoustic  wave  equation.  The  graph  below  is  presented  in 
dB-scale,  plotting  PL(dB)  versus  DEPTH  (ft). 
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Note  that  we  do  not  expect  the  solution  of  the  parabolic  equation,  multiplied 
by  an  appropriate  H^\kr),  to  agree  closely  with  the  solution  of  the  acous- 
tic wave  equation  since  the  parabolic  solution  is  just  an  approximation  of 
the  convolution.  The  introduction  of  an  additional  boundary  point  and  the 
approximation  of  a boundary  condition  by  a combination  of  normal  derivative 
and  the  parabolic  equation,  which  results  in  a second  order  ODE  of  the  ini- 
tial value  problem,  may  produce  less  accurate  results. 
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6.  CONCLUSIONS 

The  above  ODE  and  FD  methods  were  developed  in  order  to  solve  the  para- 
bolic wave  equation  with  arbitrary  bottom  and  arbitrary  boundary  conditions, 
which  the  existing  split-step  algorithm  cannot  handle.  These  methods  have 
been  shown  to  be  general  purpose  and  to  provide  the  desired  accuracy. 

At  this  stage  of  their  development,  both  ODE  and  explicit  FD  methods 
require  small  range  step-size  for  accuracy.  Implicit  FD  methods  have  favor- 
able stability,  but  the  explicit  FD  methods  do  not  have  this  property.  The 
implicit  FD  methods  are  faster  than  the  variable  step-size  ODE  methods  and 
are  equally  accurate. 

A categorization  of  the  different  environments  and  various  boundary  con- 
ditions yields  the  following  four  cases: 


CASE  1:  PLANE  PARALLEL  CONDITIONS  WITH  u = 0 at  the  BOTTOM 

SURFACE  u*0 


BOTTOM 


CASE  2:  PLANE  PARALLEL  CONDITIONS  WITH  u i 0,  ARBITRARY  AT  THE  BOTTOM 


SURFACE  u*0 
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CASE  4:  IRREGULAR  BOTTOM  WITH  NEUMANN  BOUNDARY  CONDITION 

SURFACE  oa0 
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The  applicability  of  these  methods  is  summarized  in  the  following  table, 
where  NA  stands  for  "Not  Applicable" 


^METHOD 
CASE  

ODE 

FD 

SPLIT- STEP 
(FFT) 

EXPLICIT 

IMPLICIT 

1 

X 

X 

X 

X 

2 

X 

X 

X 

NA 

3 

X 

X 

NA 

NA 

4 

X 

X 

NA 

NA 

It  is  evident  that  most  general  purpose  algorithms  used  to  solve  para- 
bolic wave  equations  are  ODE  or  explicit  FD  methods.  The  explicit  FD  methods 
are  restricted  to  the  use  of  small  step-size  in  order  to  achieve  reasonable 
accuracy;  they  are  inferior  to  ODE  methods.  In  addition,  the  built-in  cor- 
rector required  by  explicit  FD  methods  is  very  costly  in  speed  and  memory 
capacity. 
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